#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <algorithm>

#include "AggPlotter.hpp"
#include "SphData.hpp"

int main(int argc, char** argv)
{
    SphData data;
    SphData data2;
    AggPlotter plotter;
    
    plotter.readConfig("detonation__.gra");
    
    double t_start, t_end;
    
    int n_begin, n_end, n_file;
    std::string sKey;
    
    if(argc < 2){
        sKey = "q";
        return 0;
    }else{
        sKey = argv[1];
        if(argc < 3){
            n_begin = 1;
            n_end = 1;
        }else if(argc == 3){
            sKey = argv[1];
            n_begin = atoi(argv[2]);
            n_end = n_begin;
        }else{
            n_begin = atoi(argv[2]);
            n_end = atoi(argv[3]);
        }
    }
    
    char s0[100];
    int i;

    std::vector<double> xx, yy, pp, dd, vv, burn;
    std::vector<double> xx2, pp2;
    std::vector<double> xx3, yy3, pp3, dd3, vv3, burn3;
    std::vector<double> xx_3, yy_3, pp_3, dd_3, vv_3, burn_3;
    
    std::vector<double> front_xpp, front_t, front_xp;
    std::vector<double> front_xv, front_bl, front_br;
    
    std::vector<std::pair<double, int> > xxii;
    std::pair<double, int> xi;
    
    int count = 0;
    for(n_file = n_begin; n_file <= n_end; ++n_file){
        t_start = get_time();
        
        xx.clear(); yy.clear(); pp.clear(); vv.clear(); dd.clear();
        xx2.clear(); pp2.clear();
        xx3.clear(); yy3.clear(); pp3.clear(); vv3.clear(); dd3.clear(); burn3.clear();
        xx_3.clear(); yy_3.clear(); pp_3.clear(); vv_3.clear(); dd_3.clear(); burn_3.clear();
        data.clear();
        data2.clear();
        
        if(sKey == "p"){
            
            sprintf(s0, "out-%04d.dat", n_file);
            std::cout<<s0<<"\n";
            
            if(count % 5 == 0) data.dkeys.clear();
            
            data.readFile(s0);
        
            xx = data.dvalues[data.map_names["X"]];
            yy = data.dvalues[data.map_names["Y"]];
            pp = data.dvalues[data.map_names["P"]];
            dd = data.dvalues[data.map_names["D"]];
            vv = data.dvalues[data.map_names["UX"]];
            burn = data.dvalues[data.map_names["BURN"]];
            //        xx = data.getData("X");
            
            std::cout<<"Sizes: "<<xx.size()<<" "<<yy.size()<<" \n";
        
            //sprintf(s0, "agg_test-%04d.ppm", n_file);
            //plotter.plotParticles(s0, xx, yy, dd, pp);
            
            //sprintf(s0, "agg_uxtest-%04d.ppm", n_file);
            //plotter.plotParticles(s0, xx, yy, dd, pp);
        
            /* min max */
            
            std::cout<< " -- Min max --\n";
            
            int nn = xx.size();
            
            double max_x, min_x;
            min_x = *min_element(xx.begin(), xx.end());
            max_x = *max_element(xx.begin(), xx.end());
            
            double max_y = yy[0], min_y = yy[0];
            min_y = *min_element(yy.begin(), yy.end());
            max_y = *max_element(yy.begin(), yy.end());
            
            /* get band */
            
            std::cout<< " -- Passband --\n";
            
            int max_i = 0;
            
            int kFill = 0;
            if(data.dkeys.size() == 0) {
                data.dkeys.resize(nn, 0);
                kFill = 1;
            }
            
            for(i = 0; i < nn; ++i){
                if( yy[i] > data.yDown && yy[i] < data.yDown + (data.yUp - data.yDown) / 5.0 ){
                    xx_3.push_back(xx[i]);
                    yy_3.push_back(yy[i]);
                    pp_3.push_back(pp[i]);
                    dd_3.push_back(dd[i]);
                    vv_3.push_back(vv[i]);
                    burn_3.push_back(burn[i]);
                    
                    if(kFill) data.dkeys[i] = 1;
                }
            }
            
            int nn3 = xx_3.size();
            xxii.clear();
            
            for(i = 0; i < nn3; ++i){
                xi.first = xx_3[i];
                xi.second = i;
                
                xxii.push_back(xi);
            }
            
            stable_sort(xxii.begin(), xxii.end());
            
            for(i = 0; i < nn3; ++i){
                xx3.push_back(xx_3[xxii[i].second]);
                yy3.push_back(yy_3[xxii[i].second]);
                pp3.push_back(pp_3[xxii[i].second]);
                dd3.push_back(dd_3[xxii[i].second]);
                vv3.push_back(vv_3[xxii[i].second]);
                burn3.push_back(burn_3[xxii[i].second]);
            }
            
            std::cout<<"Band size: "<<xx3.size()<<"\n";
            
            //sprintf(s0, "agg_test3-%04d.ppm", n_file);
            //plotter.plotParticles(s0, xx3, yy3, dd3, pp3);
            
            /* smooth */
            //xx2 = xx;
            //pp2 = pp;
            //data.smoothData(xx2, pp2);
            /* plot */
            
            //double x1, x2, x3, x4;
            double xp, xv, xpp;
            double bl, br;
            
            std::cout<< " -- Detect front --\n";
            
            //x1 = data.detectFront(xx, pp);
            //x2 = data.detectFront(xx2, pp2);
            xp = data.detectFront(xx3, pp3);
            xpp = data.detectPeak(xx3, pp3);
            xv = data.detectFront(xx3, vv3);
            
            //sprintf(s0, "front-%04d.ppm", n_file);
            //plotter.plotLines(s0, xx, pp, x1);
            
            sprintf(s0, "re3_front3p-%04d.ppm", n_file);
            plotter.plotLines(s0, xx3, pp3, xp);
            sprintf(s0, "re3_front3pp-%04d.ppm", n_file);
            plotter.plotLines(s0, xx3, pp3, xpp);
            
            sprintf(s0, "re3_front3v-%04d.ppm", n_file);
            plotter.plotLines(s0, xx3, vv3, xv);
            //plotter.plotLines("agg_test_lines2.ppm", xx2, pp2, x2);
            
            front_t.push_back(data.tdata);
            front_xp.push_back(xp);
            front_xpp.push_back(xpp);
            front_xv.push_back(xv);
            front_bl.push_back(bl);
            front_br.push_back(br);
        
        }else{
            sprintf(s0, "out-%04d.dat", n_file);
            std::cout<<s0<<"\n";
            data.readFile(s0);
            
            sprintf(s0, "out1-%04d.dat", n_file);
            std::cout<<s0<<"\n";
            data2.readFile(s0);
            
            double err_x, err_p, err_xx, err_pp;
            
            xx = data.dvalues[data.map_names["X"]];
            pp = data.dvalues[data.map_names["P"]];
            
            xx2 = data2.dvalues[data2.map_names["X"]];
            pp2 = data2.dvalues[data2.map_names["P"]];
            
            err_xx = 0;
            err_pp = 0;
            for(i = 0; i < xx.size(); ++i){
                err_x = (abs(xx[i]) + abs(xx2[i]) > 1.0E-32)?2.0 * abs(xx[i] - xx2[i]) / (abs(xx[i]) + abs(xx2[i])):0;
                err_p = (abs(pp[i]) + abs(pp2[i]) > 1.0E-32)?2.0 * abs(pp[i] - pp2[i]) / (abs(pp[i]) + abs(pp2[i])):0;
                
                err_xx += abs(err_x);
                err_pp += abs(err_p);
            }
            err_xx /= xx.size();
            err_pp /= xx.size();
            
            std::cout<<"Errors: err_x = "<<err_xx<<" err_p = "<<err_pp<<"\n";
            
        }
        
        t_end = get_time();
        std::cout<<"<< Time step = "<<t_end - t_start<<" sec>>\n";
        
        count += 1;
    }
    
    std::vector<double> front_xp_dt, front_xpp_dt, front_xv_dt, front_bl_dt, front_br_dt;
    
    {
        data.derivate(front_t, front_xp, front_xp_dt);
        data.derivate(front_t, front_xpp, front_xpp_dt);
        data.derivate(front_t, front_xv, front_xv_dt);
    }
    
    std::ofstream ofile("re3_front.dat", std::ios_base::out);
    ofile<<"T XP XPP XV XP_dt XPP_dt XV_dt"<<"\n"; 
    for(i = 0; i < front_t.size(); ++i){
        ofile<<front_t[i]<<" ";
        ofile<<front_xp[i]<<" "<<front_xpp[i]<<" ";
        ofile<<front_xv[i]<<" ";
        //ofile<<front_bl[i]<<" "<<front_br[i]<<" ";
        ofile<<front_xp_dt[i]<<" "<<front_xpp_dt[i]<<" ";
        ofile<<front_xv_dt[i]<<"\n";
    }
    ofile.close();
    
    return 0;
}
